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Abstract 

Mechanical exfoliation of a graphite surface with an adhesive nanoasperity is studied under 
different temperatures ranging from 298 K to 2 K using classical molecular dynamics. Two 
types of the interlayer interaction are investigated. For a pairwise Lennard- Jones potential 
the complete removal of the upper graphene layer during the retraction of the nanoasperity 
occurs in the whole range of the temperatures considered. The results obtained using registry 
dependent potential, which takes into account electronic derealization contribution besides 
the van der Waals one, exhibit more pronounced temperature dependence. In this case the 
exfoliation takes place for temperatures higher than 16 K, but beginning from 8 K down 
to 2 K the system behavior manifests qualitative changes with the absence of cleavage of 
the sample. Analytical estimates combined with the results of the simulations reveal that 
the contribution of the overlap of tt orbitals of carbon atoms plays an important role in the 
exfoliation of graphite. 
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1. Introduction 

Graphite is a lamellar material which is widely used in the experiments where an atomi- 
cally flat surface is required and its fabrication is accomplished by mechanical cleavage of a 
graphite sample [U El EJ H] . The cleavage of graphite is usually considered in the mentioned 
applied context. However, understanding the detailed physics of this process and elucidating 
the influence of different factors on its behavior may be valuable both from practical and 
theoretical viewpoints. In this context it is worth mentioning that mechanical exfoliation 
was the technique which allowed the discovery of graphene [5] , a monolayer of carbon atoms 
tightly packed into a honeycomb lattice. This novel material has unusual electronic prop- 
erties and it is promising for a wide variety of applications, in particular, the creation of 
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new high-frequency electronic devices [SI Ej- In spite of the development of new methods 
for producing graphene at high yields [Bl El HO] , micromechanical cleavage or exfoliation of 
bulk graphite still remains the main technique used by most experimental groups for the 
fabrication of high-quality graphene samples [TOl ITT| H2] . Comprehensive understanding of 
this process may assist in adjusting the conditions for production of samples with the desired 
characteristics. The cleavage of graphite is also closely related to the so-called superlubricity 
which was observed during probing a graphite surface with tungsten tip of a friction force 
microscope (FFM) P, Hj. This phenomenon is characterized by a reduction of friction by 
orders of magnitude and it is attributed to the existence of a small graphite flake attached 
to the tip [H HI [131 EH]- Revealing the contributions of different factors to the formation of 
a graphitic flake, which occurs by cleavage from a graphite surface, may be valuable for the 
establishing the conditions of realization of the superlubricity phenomenon. 

In spite of practical significance of the graphite exfoliation there is a lack of its theoretical 
studies. Models of superlow friction of graphite are often based on the assumption of the 
presence of the cleaved graphitic layers [H [131 HH EE] • There are also theoretical investiga- 
tions of nanoindentation of graphite using classical molecular dynamics (MD) [161 EEO [H] or 
boundary element method [Tj5] . Diamond [TB] and virtual indentors [TTl HHJ [TO] are employed 
to probe the mechanical properties of graphite [161 [TTl EH] or to explore the formation of 
interlayer sp 3 bonds under high pressures [18J. However, repulsive interactions between the 
indentor and the sample in the works mentioned above do not allow the investigation of me- 
chanical exfoliation of graphite which could have been observed for the adhesive tips. Some 
theoretical analysis of graphite cleavage can be found in Ref. [211] where novel fabrication 
method for incorporating nanometer to micrometer scale few-layer graphene features onto 
substrates with electrostatic exfoliation is described. Numerical simulations represented in 
Ref. [20J, however, are intended only to determine the field strengths needed for performing 
the described process and do not reveal the accompanying physics. 

To fill up the gap in theoretical studies of graphite cleavage we have carried out computer 
experiments using classical MD. The considered model resembles ones described above for 
the graphite indentation, but it has the following two principal differences. The first one is 
the use of the adhesive tip. Note, that the indentation, which occurs as a consequence of 
a jump-to-contact, is very shallow and is not the target of the current investigation. The 
second difference pertains to the interlayer interaction. In the mentioned works [T6"l IT71 [T5] a 
pairwise Lennard- Jones potential (LJP) which takes into account only van der Waals (vdW) 
attraction between the layers is used. However, as studies exploiting quantum-mechanical 
techniques suggest, there is also a short-ranged electronic derealization contribution to the 
interlayer bonding of graphite [2TI [22] , and neglecting it may influence the exfoliation pro- 
cess. Thus, we performed the separate simulations using the LJ and the registry dependent 
potential (RDP) [22], which includes the mentioned orbital overlap contribution. The main 
aim of the work is to analyze the graphite cleavage under different temperatures using these 
two interlayer potentials. Temperature of the sample is one of the natural factors that 
has an impact on the interlayer cohesion in graphite, and it has been recently used in the 
solvothermal-assisted method of graphene production [9], hence indicating the need of its 
thorough exploration. The next section gives the details of the simulation setup. 

2 



2. Model 



The graphitic sample consists of three graphene layers with AB stacking (fig. [T]) which 
reflects a form of graphite. Armchair and zigzag graphene edges lie along x and y coordinate 
axes, respectively, and periodical boundary conditions are applied in the xy-p\a.ne. Each 
layer is composed of 24 x 24 honeycombs thus containing 3456 carbon atoms and the lengths 
along x and y directions are 10.082 nm and 8.731 nm, respectively. To hold the sample in 
space, the bottom graphitic layer is rigid throughout the simulations. 




Figure 1: Top (a) and perspective (b) views of the initial atomic configuration of the studied system. Green 
and cyan balls correspond to tungsten and carbon atoms respectively (all snapshots in this work are produced 
with Visual Molecular Dynamics software [25]). 

Our model is mainly approached to the experiments pertaining to the superlubricity. 
The graphite surface interacts with an infinitely hard square pyramidal nanoasperity (to 
which we also refer as the tip) which simulates the tip of FFM. The asperity consists of five 
layers of atoms parallel to the xy-plane. Particles are arranged in a perfect bcc lattice with 
constant of 0.3165 nm and this corresponds to the crystal structure of tungsten [21] . The 
tapered form is provided by adding one atomic row in x and y directions per layer when 
moving from the bottom (which is the nearest to the sample part of the asperity) to the top 
of the tip. The bottom atomic layer exposes (001) crystallographic plane and has 13 x 13 
atoms on the area a x x a y . The nanoasperity contains 1135 atoms and the total number of 
particles involved in the simulations is 11503. 

It should be noted, that the hardness of the nanoasperity may influence the exfoliation 
and, strictly speaking, for completely realistic reproduction of the experimental conditions 
the tip in the model should be able to deform. This can be achieved by exploiting one of the 
available interatomic potentials for tungsten, e. g. based on the modified embedded atom 
method |25j or its new form [26J. Nevertheless, absolutely rigid surfaces are quite often 
used in MD simulations [23 EE] and we also decided to study the system under mentioned 
approximation as the first step towards more realistic modeling. 

Nanoasperity dimensions are chosen to satisfy the fact that accordingly to the experi- 
ments the flake is assumed to attach to asperities on the tip with sizes of several nanometers 
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(see high-resolution transmission electron microscopy micrograph of the tungsten tip in 
fig. 5. 11 in Ref. [4]). The size of the nanoasperity greatly affects the exfoliation in our sim- 
ulations (this is analyzed in section [4]) and its value has been chosen to provide the most 
suitable conditions for the demonstration of the differences in system behavior with LJP 
and RDP. 

Covalent bonds between carbon atoms within the two upper dynamic graphene layers 
are described by the Brenner potential. It has the following form 



j>i 



In the current study expressions of a second-generation reactive empirical bond order (REB02) 
form of the potential [30] are used for pair-additive interactions V R (rij) and V (ry). Bond 
order function &y is chosen as in the first version of the Brenner potential (REBOl) with pa- 
rameters for potential II in Ref. [29]. The code from TREMOLO software [31] is partly used 
in calculations of cubic splines and their derivatives in the bond order term, and the inter- 
actions from Brenner potential are computed using parallel algorithm presented in Ref. [32J. 

The use of pairwise interactions from REB02 in the current model is caused by the 
fact that REBOl is incapable of proper description of any short-range hard wall repulsion 
that might prove important under high compression [301133]. However, more complex form 
of the bond order term in REB02 and hence more intensive computations forced us to 
use the b^ from REBOl because of the computational restrictions. There are, however, 
several arguments that justify its use in the context of our problem. First, let us analyze the 
roles that different contributions play in the REBO potential. The energetics of each given 
hydrocarbon structure is defined by the pairwise terms V R (rij) and V A (rij) with the latter 
modulated by the bond order function b^. The main aim of b^ is to appropriately adjust the 
energy of the atomic structure when the changes in the local atomic environment occur. This 
is accomplished by tracking the number of the nearest carbons and hydrogens and the angle 
9 between the neighboring atoms. If the atomic structure is not changed, the universal 
function b^ can be, roughly speaking, superseded by its numerical value for the current 
configuration. In our model the atomic coordination is not assumed to alter considerably, 
so we could have used the corresponding numerical value instead of b^ in eq. (11). But 
to take into account the changes in the local environment due to thermal fluctuations we 
have exploited the function b^ from REBOl, which gives the value of by close to the one 
from REB02 for graphite at the equilibrium. Substituting 6 = 120° and r = 1.42 A in 
expressions for G(0) and V B for REBOl [29] one obtains G(0) = 0.0372, % = 0.9648 and 
V B = -5.2854 eV. For REBO 2 the corresponding values are [3D]: G(cos(0)) = 0.0528, 
bij = 0.9511 and V B = —4.9861 eV, indicating the difference in binding energy between the 
two potentials only in 6 %. This is admissible for our problem because taking into account 
other approximations of the model we do not pretend to obtain accurate quantitative results. 
Note, however, that considerable deviation of an atom from the equilibrium position may 
lead to an instability as the consequence of inability of the used b^ to appropriately describe 
considerable changes in atomic coordination, resulting in the rearrangements of atoms and 
the formation of defects. As the simulations show, in our case this can be observed for 
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temperatures higher about 25 K (see section |4j). Nevertheless, the main results of our 
computer experiments pertain to lower temperatures where the instabilities do not occur 
and hence the overall conclusions of the study are not altered. 

As another arguments for the use of such a potential it should be noted that approxi- 
mations for in-plane interactions in graphite are often employed when the main concern is 
directed at the interlayer processes jU H3J HU [TJ)]. Speaking about inaccurate stiffness of 
the layers, as our estimates in section [4] show, it might have quantitative but not qualita- 
tive influence on the exfoliation and, in addition, our model provides the prospects for the 
investigation of the effect of stiffness of layers on the considered process. Lastly, a uniform 
use of just one in-plane interaction allows comparing the behavior obtained with different 
interlayer potentials, which is the aim of the present study. 

For modeling of the cleavage of graphite the crucial role may play the proper description 
of the interlayer binding [22], [33] - A pairwise Lennard- Jones potential (LJP) can describe the 
overall cohesion between graphene layers, but it is much too smooth to reflect variations in 
the relative alignment of adjacent graphitic interfaces, which is also true for other graphite 
potentials, such as AIREBO [35], that use LJ interaction for the coupling of layers. The 
reason for this is that the corrugation is mainly defined by the electronic derealization 
contribution [2T1 [22] , which is anisotropic and cannot be described in a natural way by the 
single length scale of a Lennard- Jones potential. Since during the exfoliation the upper 
graphene layer is being deformed, the anisotropic interaction may be significantly changed, 
which might greatly affect the considered process. To reveal this effect we explore two types 
of the interlayer potential. The first is a registry-dependent interlayer potential (RDP) that 
has the following form [22J : 



where A = 3.629 A" 1 , z = 3.34 A, C = 3.030 meV, A = 10.238 meV. The potential contains 
an r -6 vdW attraction and an exponentially decaying repulsion due to the interlayer wave- 
function overlap. To reflect the directionality of the overlap the function / is introduced 
which rapidly decays with the transverse distance p. The latter is defined via the distance 
Yij between pairs of atoms % and j belonging to distinct layers and the vector {k = i,j) 
which is normal to the sp 2 plane in the vicinity of atom k. In the present study is 
computed as "local" normal, i. e. as average of the three normalized cross products of the 
displacement vectors to the nearest neighbors of atom k, and this corresponds to RDP1 in 
Ref. [22] (see section [4] for more details). For long-range vdW term the cutoff distance is 
equal to r c = 2.7zo = 0.9018 nm. The presence of normals in the RDP makes it in essence a 
many-body potential which requires much more computational effort as compared to simple 
pairwise potentials. In the current study interactions only between the adjacent layers are 
considered and they are computed using our specially developed parallel algorithm based on 
linked cell lists [HUES]. 

In the second series of the simulations the interlayer energy is represented by pairwise 



Vfo.m.n,) = e-^-*)[C + f( Pij ) + f( Pji )} 
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where r is the distance between a pair of carbon atoms, values of Eqc — 2.8 meV and 
occ — 3-33 A were adjusted to obtain the interlayer binding energies and spacings close to 
ones that RDP gives, and the cutoff distance r c is the same as for the RDP. 

The tip is assumed to interact only with the upper graphitic layer. LJP with parameters 
Ewe = 0.5 eV, awe = 0.5zo corresponding to ecc, &cc m eq. (J3j) , is used for interactions be- 
tween the tungsten and carbon atoms. These values of ewe and awe provide the conditions 
for the exfoliation of the upper graphene layer under room temperature in our model [37] . 
The equations of motion are integrated using the leapfrog method [SB] with a time step 
At = 0.1 fs. 

The necessary temperature is maintained through the Berendsen thermostat coupled with 
two dynamic graphitic layers and implemented as in Ref. [31]. The thermostat accounts for 
the numerical and round errors accumulating at each integration step, and also provides the 
means to remove excess heat, which occurs as the result of work done on the system during 
approach and retraction of the tip. The implementation of the thermostat includes velocity 
rescaling by the factor /3 of the following form: 



= yi+7(^-i), (4) 

where T and T are the desired and current temperatures of the system, respectively, and 
7 G [0; 1] characterizes the rate of heat dissipation. In our simulations 7 = 0.4 is used, which 
corresponds to rather strong coupling to the heat bath in order to prevent the destruction 
of the upper graphene layer. It should be noted, that velocities are rescaled not every time 
step, but every ten time steps (or 1 fs in physical units), and in the intermediate steps the 
system is integrated without scaling. This allows to reduce the effect of the velocity scaling 
on the distribution of energy in the system. For more details see [311 138] . 

The movement of the tip proceeds as follows. After equilibration of the system during 
1 ps with the tip outside the range of interaction hung at 1.16 nm above the surface, the 
asperity was lowered towards the sample. Motion of the tip occurs by changing z-coordinates 
of tungsten atoms in increments of 0.0106 nm. The entire system is equilibrated for 40 fs in 
between displacements of the nanoasperity. After reaching the minimum distance between 
the proximal atomic layers of the two interfacing materials of about 0.108 nm, the tip is 
immediately pulled away from the surface. Mentioned quantities correspond to the rate of 
the tip movement equal to 265 m/s. The duration of the simulations is 10 or 13 ps. 



3. Results 

During 1 ps after the equilibration period at T = 298 K, when the forces between the tip 
and the sample are still zero, the average values of the interlayer distance and the binding 
energy of the upper two dynamic graphene layers are about 0.336±0.003 nm and 42±1 meV, 
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respectively. These values differ from 0.334 nm and 48 meV computed for rigid layers using 
RDP [22] by about 1 % and 15 %. The discrepancy may be attributed to the finite cutoff 
distance used in the present study, thermal fluctuations of normals and the use of local 
normals instead of semilocal ones. Nevertheless, obtained quantities are very close to the 
experimental values [22] • For the L J interaction the corresponding values at T = 298 K are 
0.337 ± 0.003 nm and 41.8 ± 0.4 meV. 

Let us consider the results obtained with RDP and LJP separately. 

3.1. Exfoliation using RDP 

Figure [2] shows force- versus-distance curves for different temperatures T obtained when 
RDP is used. These plots reflect the changes in the normal force F acting on the tip 
with a distance to the surface. In the present work this force is computed as the sum of 
z components of forces acting on tungsten atoms from the graphitic sample, and they are 
averaged over the last 10 fs of the equilibration procedure in between displacements of the 
tip. In some works additional averaging is performed on the force-displacement curves to 
filter out the noise from thermal vibrations [221 HO]- This has not been carried out in the 
present study. 

Let us consider the behavior of the model at T = 298 K in more detail (fig. [4] shows 
several snapshots of the system for this case). Following an initial slow variation of the force 
between the graphite substrate and the tungsten nanoasperity as the latter is being pulled 
toward the surface, the onset of an instability is observed, signified by a sharp increase in 
the attraction between the two (it corresponds to the negative value of the force). This is 
accompanied by a sharp decrease in the potential energy E pot of the system (fig. [3]a). The 
maximum attraction (point A in fig. [2]a and fig. |3]a, see also fig. |4]a) corresponds to a jump- 
to-contact (JC) phenomenon [381 E] which occurs via a fast process where carbon atoms 
under the asperity displace toward it in a short time span of about 0.5 ps. This phenomenon 
is further evidenced by time dependencies of the interlayer binding energy En for the upper 
two graphene layers shown in fig. [5j where it corresponds to local maximums of En observed 
in between 4 and 5 ps. JC leads to the collision of carbon atoms with absolutely rigid 
nanoasperity, which causes a sharp peak in the force-displacement curve (point B in fig. [2]a 
and fig. |3]a). Further advancement of the tip towards the sample results in an increase of 
the repulsion indicating the repulsive wall region and the indentation of the sample [38J. 
Repulsion also continues to increase during the initial stage of the retraction of the tip ( CD 
segment of the curves) which begins after 5 ps. Moderate decay of repulsion till point E is 
followed by the segment EF where the force becomes repulsive once more and it retains the 
positive sign till point G is reached. This exhibits the tendency of carbon atoms to push 
the tip upwards. Note, that the part of the force-versus-distance curve corresponding to the 
time earlier than about 8.9 ps does not unambiguously indicate the cleavage of the sample, 
as will be discussed later. A sudden change of repulsion to attraction after point G in fig. [2]a 
and fig. [3]a at about 8.5 ps is indicative of the final stage of exfoliation, where forces between 
graphene sheets at their boundaries should be overcome. The ultimate configuration has the 
completely removed upper layer (fig. ^d) corresponding to zero interlayer energy in fig. 15} 
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Figure 2: Normal force acting on the nanoasperity as it moves towards and then withdraws from the graphitic 
surface. Abscissa values correspond to the vertical distance between the rigid graphene layer and the bottom 
tungsten atomic layer. Dashed line presents the average equilibrium position of the upper carbon layer which 
is assumed to be 0.668 nm from the bottom carbon layer. 
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Figure 3: Time dependencies of the normal force acting on the tip and the potential energy of the system 
(per atom) for the maximum (a) and minimum (b) temperatures used in the simulations. 
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Figure 4: Instantaneous atomic configurations when T = 298 K corresponding to the following points in 
fig. [2]a and fig. [3)3: (a) point A, (b) D, (c) G, (d) exfoliation of the upper layer at the end of the simulation. 
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Figure 5: Time dependencies of the interlayer binding energy of the upper two graphene layers for different 
temperatures. 

The behavior of the system does not change qualitatively in the range of temperatures 
from 298 K to 16 K, as video animations show and as can be seen from fig. and fig. |5j 
Lowering the temperature leads to the decay of thermal fluctuations and it is manifested 
in the decrease of data scattering in the force-versus-distance curves. This also reduces the 
magnitudes of collisions of carbon atoms with the tip after JC and causes the reduction 
of the peak observed at point B in fig. |2]a and fig. |3]a. One can note that in the specified 
temperature range the interlayer energy Eq increases during the retraction of the tip reaching 
zero value (fig. [5]) indicating the complete exfoliation of the upper layer. 

Beginning from T = 8 K down to T = 2 K the cleavage exhibits qualitative changes. 
Let us consider the case of T = 2 K in more detail. Very small thermal fluctuations cause 
almost complete smoothing of the force dependencies on distance and time in fig. [2](i and 
fig. (fig. [6] shows several instantaneous atomic configurations for this case). When the 
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tip is lowered toward the surface the behavior is similar to the considered one for higher 
temperatures with JC manifested in a sharp minimum and the following repulsive wall 
region. Note the absence of the peak attributed to collisions of carbon atoms with the 
tip. The withdrawal part of the curves is completely different from the considered above. 
During the retraction the force acting on the tip is mainly attractive indicating the tendency 
of carbon atoms to pull the tip downwards in opposite direction to the tip movement. After 
reaching a sharp minimum (point A in fig. ^p) corresponding to the beginning of the pickup 
of carbon atoms by the tip the force remains almost constant (till point B in fig. ^p) after 
which its magnitude begins to decay reaching zero value at 11.8 ps. Time dependency of 
En in fig. [5] exhibits maximum at 8.97 ps after which it gains the initial value in contrast 
to higher temperatures where absolute value of En decreased to zero. These facts and video 
animations obtained during the simulation (see also fig. [6ji) clearly indicate that the tip 
"loses" carbon atoms, which return to the equilibrium positions of the upper layer, and the 
exfoliation does not occur. 




Figure 6: Snapshots of the system when T — 2 K corresponding to the following points in fig. |3p: (a) point 
A - 5.48 ps, (b) B - 7.5 ps, (c) C - 8.97 ps (maximum in fig. [5|, (d) D - 11.5 ps, the tip "loses" the upper 
graphene layer. 



3.2. Exfoliation using LJP 

Main results obtained when LJP is used for the coupling of graphitic interfaces are pre- 
sented in fig. [7] - 10 For temperatures above 16 K force-displacement curves (see fig. [TJa 
for T = 298 K) and time dependencies of the force and potential energies are qualitatively 
similar to the obtained with RDP (with smaller magnitude of data scattering). As video ani- 
mation sequences and fig. [8] also suggest, in the mentioned temperature range the exfoliation 
of the sample takes place. 

For temperatures beginning from 8 K down to 2 K the behavior of F and E pot during 
the time span before 8.9 ps is also qualitatively similar to the considered for RDP. For later 
times force dependencies differ qualitatively from RDP, which is manifested in the nonzero 
value of F that slowly changes in time (see fig. [9] for T = 2 K). The interlayer energy 
ultimately approaches zero value in the whole range of the temperatures (fig. fsj) , suggesting 



that the upper graphene layer is also cleaved at low temperatures (as can be seen in fig. 10 
for T = 2 K) and the tip does not "lose" the upper graphene layer 
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Figure 7: Force- versus-distance curves obtained with LJ interlayer potential for the maximum (a) and 
minimum (6) temperatures used in the simulations. 
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Figure 8: Time dependencies of the binding energy of the upper two graphene layers measured in the 
simulations employing LJP under different temperatures. 
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Figure 9: Time dependencies of the normal force acting on the tip and the potential energy of the system 
(per atom) when LJP is used for interlayer cohesion at T = 2 K. 
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Figure 10: Snapshots of the system obtained for LJP at T — 2 K corresponding to the following points in 
fig. [9j (a) point A - 8.97 ps, (6) B - 10.5 ps, (c) C - 11.5 ps, (d) D - 13 ps, the tip almost completely 
isolates the upper graphene layer. 



4. Discussion 

4-1. Qualitative analysis 

Before revealing the mechanisms leading to qualitative differences in the exfoliation ob- 
tained using the two interlayer potentials, let us more deeply analyze the potential energy 
E pot of the system and the force dependencies which also may be valuable for elucidating 
the interlayer behavior. As the binding energy En in graphite is smaller than E pot by more 
than two orders of magnitude, its contribution to the total potential energy is negligible. 
Changes in E pot are mainly defined by structural transformations in layers, by tip-sample 
interactions and by mechanical stresses occurring in the upper layer. Time dependencies of 



E pot for different temperatures obtained using RDP are summarized in fig. 11 (for LJP they 
are qualitatively similar). As was mentioned in subsection 3.1, the sharp decrease of E pot 



after 4 ps is due to JC phenomenon (see, for example, point A in fig. [2J2, fig. [3]a, fig. 12 ). It is 
followed by the fast increase, reflecting the compression of the sample, and then diminishing 
of E po t is observed till the local minimum is reached at about 5.6 ps. This corresponds to 
the initial part of retraction of the nanoasperity beginning after 5 ps. The steady increase of 
E pot in the time interval from 5.6 ps till about 8.5 ps reflects the deformation of the upper 
graphene layer during the withdrawal of the tip. Note that for temperatures higher than 
25 K (is not shown in fig. 11) the value of E pot at the moment immediately before the sharp 



increase of the potential energy at 8.5 ps is greater than its initial value before the JC. This 
is caused by the instability due to the use of the old bond order term in REBO mentioned 
in section [2] Accordingly to video animations, this results in the rearrangements of atoms in 
some parts of the upper layer into configurations different from the honeycomb one, leading 
to the formation of point defects or even of large areas with disordered structure. This 
process is more intensive under higher T and may influence the exfoliation. However, for 
temperatures lower than 25 K the honeycomb structure of the upper layer is preserved and 
this is reflected in the values of E pot . The sudden jump of E pot occurs after 8.5 ps and 
it corresponds to the moment when the interlayer binding at the boundaries of the layer 
begins to play an important role. Its onset does not depend on T indicating that this jump 
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Figure 11: Time dependencies of the potential energy of the system for different temperatures obtained 
using RDP. 



is defined by the geometry of the model and the tip-sample interaction, but not by the 
potential functions used for the interactions between carbon atoms. As video animations 
show, at this moment the tip abruptly "loses" some part of carbon atoms which results 
in sharp increase in E votl and the corresponding force dependencies show the decrease of 
magnitude of F (see, for example, point G in fig. [3]a, jump in fig. [3]£> before point C, and 
point B in fig. 12). This is a crucial moment for exfoliation and it is important that in our 



model this moment occurs immediately before the retraction stops at 8.96 ps. If the tip had 
been withdrawn further, it would have "lost" atoms also for LJP at lower temperatures. 

The above consideration indicates that force-displacement curves do not unambiguously 
reflect the possibility of the exfoliation for times less than 8.96 ps, as at this period they 
mainly reflect the state of carbon atoms under the tip but not of the whole layer. This 
explains the qualitatively similar behavior of force curves for the two potentials (see also 



fig. 12). For times greater than 8.96 ps cleavage corresponds to zero force (when no atoms 



interact with the tip, point C in fig. 12) which is observed for RDP. Comparing the depen- 
dencies of -E^t for T = 16 K and T = 2 K in fig. 11 and corresponding curves for F in fig. |2]c 
and fig. [2}i, respectively, may appear to be inconsistent at a first glance, as the force curves 
are very different in contrast to E pot . However, accordingly to the above discussion, the 
discrepancy may be ascribed to the different behavior of carbon atoms under the tip during 
retraction. At the higher temperature they tend to move upwards after each withdrawal 
step, thus pushing the tip in this direction, but for the lower temperature due to smaller 
mobility they pull the asperity in the opposite direction. It should be noted that decrease of 
E pot with diminishing the temperature may be explained using the technique applied below 



in subsection 4.2 for interlayer energy, but it is out the scope of the current investigation. 
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Figure 12: Force-versus-distance curves obtained using LJP and RDP at T = 8 K. 



Time dependencies of the interlayer binding energy En of the upper two graphene layers 
(fig. [5] and fig. [8]) provide the means to unambiguously treat the interlayer processes. The 
apparent interaction of the upper layer with the tip is manifested in the maximum after 
about 4 ps due to JC, although small but nonzero value of the force occurs after about 
2 ps. As fig. [5] and fig. [8] suggest, En diminishes with the decrease of T. Averaging of 
En and of the interlayer distance d during the time span from 1 ps (after equilibration) 
till 2 ps (when the tip does not interact with the sample) shows that for LJP mentioned 
tendency is observed in the whole temperature range, suggesting that thermal expansion of 
the interlayer spacing occurs with increase of T. For RDP this behavior is also observed 
for temperatures higher than 16 K. However, En and d do not decrease for RDP when T 
diminishes from 8 K to 2 K and even the opposite trend is observed. This may be ascribed 
to a rather small averaging time interval not enough to obtain the true average values, as 
for RDP much higher fluctuations of En are observed. On the other hand, these results may 
point out that although at high temperatures thermal expansion greatly facilitates cleavage, 
at low T its significance diminishes and another contribution plays a crucial role. 

The withdrawal parts of curves for En at high temperatures are monotonic and quali- 
tatively similar. They may suggest that defects in the upper layer and thermal expansion 
facile the cleavage process. For LJP with decrease of T the slope of the curves diminishes, 
indicating the slowing of the exfoliation, but the form of the dependency is preserved down 
to 2 K (see also fig. 13). In contrast, for RDP lowering T leads to qualitative changes. At 



about 50 K the curve becomes nonmonotonic at about 8.5 ps, which corresponds to a mo- 
ment of a sharp jump in E pot (see the above discussion). Further decrease of T transforms 
bending of the curve into a plateau and ultimately into a maximum. These results clearly 
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Figure 13: Time dependencies of the binding energy of the upper two graphene layers obtained using RDP 
(dashed lines) and LJP (solid lines) at low temperatures. 

suggest that for RDP a potential barrier exists, which should be overcome at the moment 
after about 8.5 ps to provide cleavage of the sample. The barrier is not apparent at high 
temperatures but manifests itself with decrease of T. In our model, beginning from 8 K 
the magnitudes of the tip-sample interaction and of the stresses in the upper layer are not 
enough to overcome the barrier and thus the exfoliation does not occur. 



4-2. Phenomenology 

Now we make some analytical estimates which may help to conceive the appearance of 
the energy barrier for RDP and to explain the observed differences for the two potentials. 
Note, however, that we do not pretend to obtain accurate quantitative results and the main 
aim is to reveal the main trends. 

At first, let us analyze the behavior of the RDP under different temperatures. The main 
feature of this potential is function / (see eq. ((21) ) of the following form [22J: 



f(p) 



Cn + Co 



+ c 4 ( 



exp 



-V 

SJ 



(5) 



Here Co = 15.71, Ci = 12.29, C4 = 4.933 are measured in meV and S = 0.578 A. Quantity / 
reflects the directionality of the overlap of 7r orbitals and makes the dominant contribution 



to the repulsive part of the RDP. It rapidly decays with the transverse distance p (fig. 14) 
which is defined using local normal in the vicinity of an atom k: 



Pi 



r 2 



n, ; r, 



)V 



J I 



r 2 



(n 



(6) 
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For a form of graphite atoms in distinct layers can appear one under another, or an 
atom in one layer can be located under the center of a honeycomb of another layer. We 
consider only the latter situation, for the former the trends should be qualitatively similar. 
It is sufficient to study the contributions to / only of the first neighbors as the transverse 
distance p to the second and farther neighbors is more than 2 A, and thus their contribution 
to / is negligible (see fig. [lift). For static layers the distance between the mentioned two 
atoms from distinct layers is r = a/ ' Zq + a 2 , where a = 1.42 A is the interatomic distance 



in a graphite layer (see inset in fig. 14b). In dynamic layers radius- vector connecting such 
two atoms can be written as follows: 

r(i) = r„ + u(*), (7) 

where u(t) is a stochastic quantity that changes in time due to thermal fluctuations. Ther- 
mal motion of atoms also leads to fluctuations of normals, that are used to determine the 
orientation of tc orbitals. Therefore, the angle 9 between the vectors r and n (fig. [i~4|a) also 
fluctuates and we can write 

e{t) = Q + e{t), (8) 

here 9$ — arcsin(a/ro) corresponds to static layers and 9(t) changes in a stochastic manner. 
Substituting (17]), ^ in ([6]), and applying averaging by sufficiently long time interval gives 
the following expression for a mean square of p (we omit indices for simplicity): 

(p 2 ) = ( r 2 ) - (r 2 cos 2 9) = (r 2 + u 2 + 2r u) 

(r 2 + u 2 + 2r u) cos 2 (9 + 9)^. (9) 
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Since the normal of an atom is defined using the nearest neighbors in the layer, and vector 
u pertains to atoms in different layers, these quantities can be assumed to fluctuate inde- 
pendently. Thus, averaging the expressions containing u and 9 can be performed separately. 
With this in mind, applying standard trigonometric identity and taking into account a small 
magnitude of 9 and that (i"ou) = 0, one obtains 



(p 2 ) = r 2 + (u 2 ) - r 2 ( (cos 2 9 - 9 sin 29 + 9 2 sin 2 9 
- (u 2 (cos 2 9 - ^sin 29 + 9 2 sin 2 9 



(10) 



Employing the equality (9) = we have: 

<p 2 >=(;r 2 + < M 2 >-r 2 (^ 2 )-< M 2 >(^ 2 ))sin 2 ^. 



11^ 



From the thermodynamic relation C (u 2 ) /2 = 3ksT '/2 follows that (u 2 ) = aT, where a = 
3ks/C, k-Q is the Boltzmann constant and C is an effective spring constant. Substituting 
values of carbon atomic mass M = 1.994 • 10~ 26 kg and of frequency u> = 10 14 rad/s in 
C = Moo 2 , one obtains C = 199.4 N/m and a = 2.1 • 10~ 5 A 2 /K. Analogously assuming 



that (9 2 ) = 0T and \/ (9 2 ) = 15° at 300 K, we have (3 = 2.25 • 10" 4 rad 2 /K. Eq. (|ll|) takes 
the form 



(p 2 ) = a 2 + (aT - f3r 2 T - a(3T 2 )a 2 /r 2 . 



(12) 



Inset in fig. 15 plots the dependency a/ (p 2 )(T). As can be seen, it decays with T. 
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Further step is to express analytically (/) through (p 2 ). One way to do this is to expand 
exp [— (p/5) 2 ] into series by powers of p/5 near zero point and then substitute expression 



from eq. (12). However, since p/5 is not small this leads to alternating series with increasing 
terms, and therefore we cannot terminate the series at finite number of terms. Direct 
integration for averaging also does not allow expressing analytically (/) through (p 2 ) and 
the numerical techniques should be involved. To avoid this, we simply substitute (p 2 ) in 
eq. (p)| to obtain very rough estimate of the dependence (f)(T). The result is shown in 
fig. [if 

For vdW interaction neighbors farther than the first one should be taken into account to 
obtain accurate energy values. This is the reason for relatively large cutoff distance used in 
the simulations which covers about 6 or 7 neighbors. Nevertheless, the major contribution 
makes the first neighbor and for estimative aims it is sufficient to consider it. It can be 
shown that with accuracy to moments higher than the second the expression (r 6 ) ph (r 2 ) 3 
is true. Using eq. (fin), we obtain for the temperature dependence of magnitude of vdW 
interaction from eq. (pi) : 

4-6 A ~6 A r 6 

iv \ - ~ - h *\ 

[V vdw) - + u)6) ~ ^ + u2)3 - { ^ + aT)r [16) 



The results are summarized in fig. [15} As can be seen, / increases with temperature, 
that corresponds to greater repulsion due to the interlayer wave-function overlap, and vdW 
attraction decreases, indicating thermal expansion of the interlayer separation. It can be 
noted that the growth rate df/dT is by about an order of magnitude larger than the rate 
of vdW decaying with T. 

Let us analyze the formation of the potential barrier. It is defined by the sizes of the 
contact surface of the tip a x and a y , the dimensions of the sample L x , L y , the ultimate tip 
height and the rate of its retraction v, by stiffness of the upper layer and by the energy of the 
tip-sample interactions ewe- We consider the situation which is observed in the simulations 
for low temperatures, when for the given ewe the rate v has such a value that the upper 
layer is neither completely cleaved nor is "lost" when the tip reaches the height h (above the 



equilibrium position of the upper layer, see fig. 16) at about 8.6 ps corresponding to jumps 
in E pot immediately before the tip stops. Three regions with different contributions to the 



interlayer interaction can be marked out (fig. 16) 



1. Region I - atoms are not effected by the tip. Here both contributions - repulsive from 
the Ti orbital overlap V* and attractive vdW Vj dW are presented (these quantities stand 
for average values of the interlayer energy). 

2. Region II - both contributions V™ and V^l w exist but the magnitude of V^ 1 decreases 
to when approaching to region III. 

3. Region III - only vdW attraction V^l w is presented. 

For the fixed values of v, ewe, h and of the interlayer potential the dimensions of these 
zones are defined by the geometry of the system and the stiffness of the layer. We denote 
by i?2 the distance from the center of symmetry of the layers to the boundary between the 
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Figure 16: Sectional view of the system when T = 8 K at 8.6 ps illustrating the three regions with different 
contributions to the interlayer binding. 



regions I and II, and by R3 the one to the boundary between regions II and III (boundaries 
are assumed to be the circles on the surface of the layer). From fig. 16 we can conclude that 

R 2 = a x /2 + htamp, (14) 

here we use the fact a x = a y in our model. The quantity R3 corresponds to some critical 
distance r cr between the layers, where p = p cr = 2 A and therefore / ~ 0. As the angle 
between the normal to the upper layer and r cr is tt/2 — ip, hence p 2 = r 2 r — r 2 cr cos 2 {n /2 — ip) 
and ultimately r cr = p cr / cosy? ~ 4.7 A for p> = 65° (the value chosen only for estimative 
purpose). From fig. 16 one can note, that 

R 3 = R 2 - (r cr - z ) tan ip. (15) 

Having defined the dimensions of the regions we can calculate the interlayer energy for 
the considered time moment. Denoting the number of atoms per unit surface area in the 
layer by n and assuming the absence of defects (which is the case for low temperatures in 
the simulations), the numbers of atoms Ni, Nu Nui in each region are defined as follows 

R2 _ R2 

A/ = n(L x L v - nRl), N H = nir^-^ — -, 

sin <p> 

N IU ~ n{a x a y + tt 6 , 2 -). 16 

sin ip 

For the binding energy of the upper two layers we have the following expression 

E a = Nj(Vj + Vj dW ) + N n (V" + V» w ) + N m V v %. (17) 



Eq. (17) presents a wrapped form of the approximation for the potential barrier that should 
be overcome to provide the cleavage of the layer. We analyze this expression qualitatively. As 
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was shown above, under the lowest temperature the repulsive contributions from the orbital 
overlap , V^ 1 have the smallest possible value, and vdW attraction between the layers 
is the largest. Therefore, in the competition of the three types of interactions: vdW tip- 
sample interaction and tt orbital overlap that tend to separate the layers and vdW attraction 
between the layers, wins the latter and the cleavage does not take place. This situation for 
the chosen parameters of our model remains up to 8 K, although the interlayer barrier 
considerably diminishes (fig. [5]). Beginning from 16 K due to fast growth of / with T the 
two interactions tending to isolate the layer begin to prevail over vdW interlayer attraction 
(which has been diminished due to thermal expansion), the magnitude of the barrier reaches 
the value that can be overcome and the exfoliation occurs. The above discussion shows that 
it is this anisotropic orbital overlap contribution that plays an important role in the behavior 
of the model due to its fast growth with T. The pairwise interaction gives minor contribution 
as the result of the weaker temperature dependence. 

It should be noted that another combination of the parameters of the model can result in 
different scenario from the observed above. For example, for much larger sample the cleavage 
could not have been occurred, since the first region would have much larger size. Stiffness of 
the layers also plays an important role as it defines the angle ip (when the other parameters 
are fixed) and hence the dimensions of the marked regions. Lower stiffness leads to smaller 
ip and enlarges the first zone, thus worsen the conditions for cleavage. The use of tips with 
reduced a x and a y will also lead to decreased the third and the second regions and hence to 
the probable absence of exfoliation. This firmly suggests that in simulations where another 
in-plane potential is employed or in experiments the observed processes would take place 
in distinct temperature range. There also might be the need to adjust other parameters to 
obtain the results described in the current work. 

5. Conclusion 

Computer experiments presented in the paper reveal the influence of the temperature 
on the exfoliation of a graphitic sample for two different types of the interlayer potential. 
The main result of the simulations is that the inclusion of the anisotropic contribution 
accounting for tt orbital overlap into the interlayer potential can qualitatively change the 
kinetics of cleavage under low temperatures, although the potentials can give almost the 
same interlayer cohesive energy. In our model this is manifested in the absence of the 
exfoliation below 8 K when RDP is used, while for LJP the upper graphene layer is isolated 
in the whole range of the considered temperatures. Some analytical estimates have been 
carried out that qualitatively explain the behavior observed in the simulations and define the 
contributions of different geometrical parameters of the system to the considered processes. 
The results obtained indicate the need of the experimental verification of the role of orbital 
overlap in the interlayer cohesion of graphite and the correctness of its description by the 
RDP. Our model provides a sketch for the possible experimental setup, which can be based 
on the electrostatic exfoliation technique presented in Ref. [20J as it allows adjusting the 
magnitude of the tip-sample interactions. 
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